1. Field of the Invention
This invention relates to seismic surveying, image processing, and more particularly, to illumination analysis of subsurface seismic horizons.
2. Description of the Related Art
Generally, marine seismic surveys are conducted by towing an energy source (emitter) and seismic receivers (detectors) behind a vessel. The source emits an acoustic wave in the water that travels downward to the underlying elastic earth. The sound wave is reflected by the various earth formations, and travels upward to the receivers. The seismic receivers (detectors) convert the acoustic signals, such as pressures, to electrical signals and those signals are thereafter recorded in various formats. Computers are used to analyze these data on board the vessel in real time and/or at on-shore office locations. The subsurface earth formations are reconstructed from this analysis, usually through images with color applied in the course of the analysis in order to differentiate the earth formations and strata. These techniques are used by geophysicists and geologists to locate particular earth formations in order to better predict sites for successful exploration wells and/or to better develop existing fields and wells.
A significant amount of exploratory work has been done in the seismic field, both onshore and offshore. The significance of finding accurate predictive tools has been growing as the exploration and production costs have continued to grow. The costs of production are prohibitive as operators drill in deeper waters and more hostile environments. Similarly, offshore seismic work has also become more expensive and difficult as this work is performed in deeper waters and more hostile environments. Currently the state of the art is represented by the following:
Assuming a spike source i.e. a source with a flat power spectrum, then the recorded data D for a wave that is emitted from the source S, propagated through the water and down into the earth, reflected at reflector r, propagated back to the surface, and recorded at the receiver G is represented by equation (1)
D(Sxe2x88x92 greater than G)=W+(Sxe2x88x92 greater than r)**R(r)**Wxe2x88x92(rxe2x88x92 greater than G) xe2x80x83xe2x80x83Eq. (1) 
(** denotes convolution operation. W+ is a forward propagating wave function that moves the seismic energy from the source to the specular reflection point r of the reflecting surface. R(r) is the reflectivity at location r. Wxe2x88x92 is a backward propagating wave function.)
F(r)=W+(Sxe2x88x92 greater than r)**W(rxe2x88x92 greater than G) xe2x80x83xe2x80x83Eq. (2) 
F(r) represents the combined propagation effects for the downward and the upward waves that reach the subsurface point r and is commonly referred to as the overburden effect between the source and receivers.
Function D(Sxe2x88x92 greater than G) is the data received by receivers. W+(Sxe2x88x92 greater than r) and Wxe2x88x92(rxe2x88x92 greater than G) are wave functions that depend on the underlying earth model. This earth model consists of a topological description of the subsurface geological formations, i.e. their relative locations, shapes and orientations, and the sound speed within these formations. R(r) is a function related to the contrast in intrinsic properties of the subsurface formations on either side of the reflecting surface. Function F(r) depends on the geometrical relationship between the source, the subsurface volume and receivers as well as the sound speed distribution in the subsurface volume. From the data function D and the overburden function F, one can determine function R which is the goal of seismic surveys.
To solve the wave equation by computing the functions W+(Sxe2x88x92 greater than r) and Wxe2x88x92(rxe2x88x92 greater than G) to derive function R(r) requires very intensive and costly computations. A highly simplified approach for computing a rough approximation of W+(Sxe2x88x92 greater than r) and Wxe2x88x92(rxe2x88x92 greater than G) is to use ray-tracing. The resulting function R(r), a function of a 3-D volume is commonly presented as a set of vertical or/and horizontal images, thus a 3-D object can be presented and analyzed more conveniently in two dimensions.
While deriving R(r) is the goal of any seismic experiment, the function F(r) is also of significant importance and use. When F(r) is computed by adding the contributions from all source and receiver pairs in a seismic survey then it represents the extent to which the subsurface has been illuminated given the particular seismic geometry that is being surveyed. The illumination information contained in F(r) helps to explain why a portion of the subsurface seismic image R(r) may be weak, noisy, unclear, etc. The most commonly used approach for estimating subsurface illumination involves the use of two-point ray tracing. Typically, a bundle of rays are shot from each surface shot location, propagated through a predefined earth model, impinge upon the horizon (a subsurface) which is the subject of illumination analysis, and reflect back to the surface where they are accepted or rejected based on the proximity of their emergence point to predefined receiver locations. The illumination maps generated by this method are either a direct or a weighted sum of the reflection points (xe2x80x9chitsxe2x80x9d) per unit surface element of the reflecting surface.
While the above approach is conceptually and computationally straightforward, in reality it has a number of shortcomings. These shortcomings can result in costly operations that do not produce their intended results.
Ray tracing is very sensitive to the high spatial frequencies of the velocity model as is the case for steep, high velocity contrast boundaries. For example, in the presence of salt, small but discontinuous changes in the salt sediment interface geometry (i.e. slope) quite often result in xe2x80x9cemptyxe2x80x9d receivers. In such cases, iterative shooting of infill rays is being used to ensure that rays are captured at an adequate number at receiver locations for each shot. However, frequently, infill ray shooting fails to converge, and as a result, a large number of receivers have no rays associated with them.
To overcome the problem of lack of xe2x80x9ccapturedxe2x80x9d rays, one is forced to smooth the model using mathematical algorithms, such as Boxcar smoothing algorithm, occasionally at unacceptable levels.
The ray tracing approach is generally limited to a single arrival per receiver, although in the presence of complicated subsurface bodies such as salt, this is not generally the case.
Using one arrival per receiver implies that reflection points for successive receivers may jump around as the corresponding rays may belong to different branches of an expanding wavefront.
Using single reflecting points as a way of producing illumination maps is in general a very poor approximation. In reality, the energy that is constructively interfering to produce a reflection event from an interface is reflected from a rather large area, the so called Fresnel zone. The Fresnel zone is, in addition to the dominant frequency, dependent upon the 3-D geometrical relationship between the impinging wavefront curvature and the local reflecting surface curvature. As a result, it is spatially highly variable and can only be poorly approximated by using boxcar smoothing techniques of the reflection point xe2x80x9chitxe2x80x9d map.
xe2x80x9cHitxe2x80x9d maps by themselves are of value only for those cases where the overburden effects on ray amplitudes (transmission losses, geometrical divergence) are more or less constant or very slowly varying and therefore can be ignored. In the case of subsalt illumination, xe2x80x9chitxe2x80x9d maps tend to be meaningless unless these overburden effects are taken into account. Computing such effects from ray tracing is in theory feasible but in practice extremely hard to accomplish with any degree of confidence.
The compounded effect of all the factors outlined above limits the utility of the illumination analysis, as it is currently practiced, to very simple models. It is desirable to have a mathematical modeling method that improves the utility of the illumination analysis.
The present invention employs a wave equation-raytracing hybrid method to generate an illumination image. In one embodiment of the present invention, the rays are one-way traced from surface to subsurface; the reflection points are found by using Fermat""s principle; and the amplitudes are calculated by wave function.